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w/'of  SPECIAL 


I . INTRODUCTION 


We  shall  be  concerned  with  iteratively  determining  the  N-vector  u 
of  a large,  sparse  linear  system 


where  A is  a real,  symmetric,  positive  definite  matrix  of  order  N and  b 
is  a given  N-vector.  Such  systems  arise  in  the  finite  difference 
solution  of  elliptic  boundary  value  problems.  Particularly,  we  shall 
develop  an  adaptive  scheme  based  on  the  symmetric  SOR  (SSOR)  iterative 
method  with  Chebyshev  acceleration.  Related  work  which  has  recently 
appeared  includes  Axelsson*,  Hayes  and  Young, ^ and  Young. 3»4,5,6 


II.  BASIC  METHOD 


By  defining 


where  D = diag  (A)  and  L and  U are  strictly  lower  and  upper  triangular 
matrices,  respectively,  we  may  replace  the  system  (1)  by 

u = Bu  + c. 


Axelsson,  0.  (1972),  "A  Generalized  SSOR  Method,"  BIT  13,  443-467. 

Hayes,  L.  J.  and  D.  M.  Young  (1977),  "The  Accelerated  SSOR  Method  for 
Solving  Large  Linear  Systems:  Preliminary  Report,"  Report  CNA-123, 
Center  for  Numerical  Analysis,  The  University  of  Texas,  Austin,  Texas 

Young,  D.  M.  (1974a),  "On  the  Accelerated  SSOR  Method  for  Solving 
Elliptic  Boundary  Value  Problems,"  in  Lecture  Notes  in  Mathematics, 

A.  Dold  and  B.  Eckmann  (eds.),  Vol.  363,  Conference  on  the  Numerical 
Solution  of  Differential  Equations,  Dundee- L973,  ftf.  A . Watson,  ed.), 
Springer-Verlag,  New  York,  195-206. 

Young,  D.  M.  (1974b),  "Solution  of  Linear  Systems  of  Equations," 
in  Numerical  Solutions  of  Partial  Differential  Equations,  J.G. 

Gram  (ed.) , D.  Reidel  Pub.  Col.,  Holland,  35-54  (Proceedings  of 
Conference  "Advanced  Study  Institute  on  Numerical  Solution  of  Partial 
Differential  Equations,"  Kjeller,  Norway,  August  20-24,  1973). 

Young,  D.  M.  (1974c),  "Stopping  Criteria  and  Adaptive  Parameter 
Determination  for  Iterative  Methods  for  Solving  Large  Linear  Systems,' 
Proceedings  of  the  Gatlinburg  VI  Symposium  on  Numerical  Algebra, 
Munchen,  Germany,  December  15-22,  1974. 

Young,  D.  M.  (1977),  "On  the  Accelerated  SSOR  Method  for  Solving 
Large  Linear  Systems,"  Advances  in  Mathematics  23,  215-271. 


The  SSOR  method  (Sheldon)7  is  then  defined  by  forming  a single  SSOR 
iteration  from  a forward  SOR  iteration  followed  by  a backward  SOR 
iteration;  that  is,  for  n = 0,  1,  2,  ...  we  set 


where  ul  7 is  an  arbitrary  initial  approximation  to  the  solution  u,  and  w 
is  a real  relaxation  parameter  such  that  0<co<2. 


Elimination  of  u 


where 


Note  that  L ^ corresponds  to  the  familiar  SOR  iteration  matrix.  The 
backward  SOR  operator  U ^ is  defined  analogously. 

If  storage  for  an  extra  N-vector  is  provided,  the  work  required  for 
one  SSOR  iteration  may  be  reduced  to  about  the  work  necessary  for  a single 
SOR  iteration.  The  work-saving  technique  is  due  to  Niethammei8  and 


Sheldon,  J.  W.  (1955),  "On  the  numerical  Solution  of  Elliptic 
Difference  Equations,"  Math.  Tables  Aids  Comput.  9,  101-112. 


Relaxation  bei  Komplexen  Matrizen 


Niethammer,  W.  (1964) 
Zeitsch.  86,  34-40. 


is  described  in  Benokraitis 
by  Conrad  and  Wallach.11 


1 


and  Young. ^ The  method  has  been  rediscovered 

-S’ 

The  SSOR  method  converges  if  S(S^),  the  spectral  radius  of  the  iteration 

matrix  S^,  is  less  than  1,  which  holds  if  0<w<2  and  A is  positive  definite. 

The  rate  of  convergence  is  governed  by  the  ordering  of  the  equations  and  by 
the  parameter  u.  Assuming  the  natural  ordering.  Young*''*'6  has  shown 
that  for  a certain  discrete  generalized  Dirichlet  problem  one  can 
choose  a "good"  u>  depending  on  bounds  for  the  eigenvalues  of  B and  LU  so 
that  the  SSOR  method  converges  with  the  same  order-of-magnitude  as  the  SOR 
method.  For  a finite  difference  discretization  with  mesh  size  h,  the 
number  of  iterations  required  for  both  methods  increases  like  h”  . 


Therefore,  even  by  employing  Niethammer's  work-saving  scheme,  there 
is  little  justification  for  using  SSOR.  However,  the  SSOR  method  can  be 
accelerated  by  an  order-of-magnitude  by  means  of  Chebyshev  semi-iteration 
since  the  eigenvalues  of  the  matrix  S are  real  and  nonnegative.  (Chebyshev 

u n 13 

semi-iteration  was  first  studied  by  Varga  " and  Golub  and  Varga.) 

This  approach  is  precluded  for  SOR  with  optimum  to  = o>  since  many  of  the 

l?  h 

eigenvalues  of  L are  complex.  (See  Varga  “ and  Young  .)  Also, 

there  is  no  improvement  when  semi-iteration  is  applied  to  SOR  with  1<io<(o.  . 

15  k 

(See  Kincaid  .)  For  accelerating  the  Gauss-Seidel  method  (SOR  with 

to  = 1),  see  Sheldon*®  and  Young.*1* 


9Benokraitis,  V.J.  (1974),  "On  the  Adaptive  Acceleration  of  Symmetric 
Successive  Overrelaxation,"  Doctoral  Thesis,  University  of  Texas,  Austin. 


10Benokraitis,  V.J.  (1976),  "An  Improved  Iterative  Method  for  Optimizing 
Symmetric  Successive  Overrelaxation,"  in  ARO  Report  76-3,  Proceedings 
of  the  1976  Army  Numerical  Analysis  and  Computers  Conference,  133-140. 

^Conrad,  V.  and  Y.  Wallach  (1977),  "A  Faster  SSOR  Algorithm,"  Numer.  Math. 

27,  371-372. 

1 2 

Varga,  R.S.  (1957),  "A  Comparison  of  the  Successive  Overrelaxation  Method 
and  Semi- Iterative  Methods  Using  Chebyshev  Polynominals,"  J.SIAM,  5,  39-46. 

13Golub,  G.H.  and  R.S.  Varga  (1961),  "Chebyshev  Semi-iterative  Methods, 
Successive  Overrelaxation  Iterative  Methods,  and  Second-order  Richardson 
Iterative  Methods,"  Numer . Math. , Parts  I and  II,  3,  147-168. 

14 

Young,  D.M.  (1971),  Iterative  Solution  of  Large  Linear  Systems,  Academic 
Press,  New  York. 

*5Kincaid,  D.R.  (1974),  "On  Complex  Second-degree  Iterative  Methods," 

SIAM  J.  Numer.  Anal . 11,  211-218. 

16Sheldon,  J.W.  (1959),  "On  the  Spectral  Norms  of  Several  Iterative 
Processes,"  J.  Assoc . Comput . Mach.  5,  39-46. 
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III.  ACCELERATED  METHOD 


The  optimum  serai-iterative  method  based  on  SSOR,  denoted  by  SSOR 
SI,  is  defined  by 


where 


In  order  to  apply  the  SSOR-SI  method,  we  must  determine  the  two  parameters 
w and  SfS^).  Some  a priori  methods  for  obtaining  these  parameters  are 
discussed  by  Habetler  and  Wachspress,  Evans  and  Forrington,*8  Young3,4,6 
and  Benokraitis.9**0  Since  finding  these  parameters  may  take  as  much  work 
as  solving  the  original  problem,  we  are  led  to  consider  adaptive  tech- 
niques which  approximate  the  parameters  and  at  the  same  time  improve  the 
solution  of  the  linear  system. 


IV.  FOUNDATION  FOR  ADAPTIVE  METHOD 


We  begin  by  characterizing  the  eigenvalues  of  Sco  in  terms  of  certain 
inner  products.  This  result  is  due  to  Habetler  and  Wachspress.*7  See 
also  Young. 

THEOREM  1 . Let  X be  an  eigenvalue  of  Sw  where  0<w<2  and  let  v be  an 
associated  eigenvector.  Then  X may  be  represented  by 


Habetler,  G.  J.,  and  E.  L.  Wachspress  (1961),  "Symmetric  Successive 
Overrelaxation  in  Solving  Diffusion  Difference  Equations,"  Math. 

Comp.  15,  356-362. 

‘Evans,  D.  J.,  and  C.  V.  D.  Forrington  (1963),  "An  Iterative  Method  for 
Optimizing  Symmetric  Successive  Overrelaxation,"  Comput . J.  6,  271-273 


(4) 


where 


X = 1 - io(2-u>) 


l-o 


1-uia+w  6 


= <Kw,v) 


(5) 


(v,DBv) 

(v,Dv) 


= 1 


v.DLUv) 


THEOREM  2.  The  representation  <K<o,v)  given  by  (4)  for  any  vector  v M is 
a Rayleigh  quotient  with  respect  to  the  vector 


w = (I 


-0)11)0^ 


and  the  symmetric  matrix 


where 


3u  = (I-wU^Sjf^fl-uU)"1 


U = D^UD-*5  . 


That  is. 


(w,S  w) 

v U) 


= (w,w)“  • 


Furthermore,  S is  similar  to  S and 

U)  0) 


(6) 


♦ (»,v)  <S(3u)-S(SJ. 


Proof:  See  Benokraitis  (1974). 

We  emphasize  that  (6)  holds  for  any  nonzero  vector  v,  not  just  for  eigen- 
vectors of  S . However,  the  closer  we  approach  a fundamental  eigenvector, 
the  closer  w we  shall  be  able  to  determine  S(S  ) from  <Ho>,v)  given  by  (4). 
Therefore,  it  would  be  fortunate  if  somehow  we  could  determine  the  fundamental 
eigenvector  without  deviating  from  the  path  of  improving  the  approximate 
solution  of  (1).  A clue  leading  to  the  desired  situation  is  contained  in  the 
following  theorem  (cf.  Young  (1974c)). 5 


9 


— 


THEOREM  3,  The  psudo-residual  vector 

(7)  6(n)  = S u(ni  + k - u(n) 

b)  b> 

where  u^  is  the  latest  SSOR-SI  iterate,  satisfies 

(8)  6(n)  = P (S  )£(0)  . 

n or 


Here 


P (S  ) = 

n a) 


r-(4> _I) 


th 


where  T^(x)  is  the  n degree  Chebyshev  polynomial  defined  by  the  three-term 
recurrence  relation 


Tq(x)  = 1,  Tj (x)  = x 

Tn+1  = 2xTn^x)  ' Tn_i  (x) » n >_  1 . 

o 

Proof:  See  Benokraitis  (1974). 

We  note  that  if  Pn(^w)  is  replaced  by  S^n  then  (8)  reminds  us  of  the  power 

method  for  computing  the  dominant  eigenvector.  With  this  motivation,  the  next 
theorem  comes  as  no  surprise. 

THEOREM  4.  The  pseudo-residual  vector  6*-0-1  given  by  (7)  converges  in  direction 
to  the  eigenvector  associated  with  the  eigenvalue  S(S  ) as  n tends  to  infinity. 

b>  J 

Proof:  See  Benokraitis  (1974) . 9 Also  compare  with  Diamond  (1971)19  and 
Hageman  (1972).  0 


19Diamond,  M.A.  (1971),  "An  Economical  Algorithm  for  the  Solution  of  Finite 
Difference  Equations,"  Report  UIUC  DCS-R-71-492,  Department  of  Computer 
Science,  University  of  Illinois  at  Urbana-Champaign,  Urbana,  Illinois. 

20Hageman,  L.A.  (1972),  "The  Estimation  of  Acceleration  Parameters  for  the 
Chebyshev  Polynomial  and  the  Successive  Overrelaxation  Iteration  Methods," 
Report  WAPD-TM-1038,  Bettis  Atomic  Power  Laboratory,  Westinghouse 
Electric  Corp.,  Pittsburgh,  Pennsylvania. 
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. i ii  ■ mi  f '..xi. 


By  Theorem  4 it  is  possible,  then,  to  determine  the  fundamental  eigen- 
vector with  a little  additional  effort  by  computing  the  pseudo-residual  vector. 
However,  since 


jC")  . s u00  . k . utn) 


' (n)  (n) 

= u - u 


one  SSOR 
iteration 


latest  SSOR-SI 
iteration 


and  since 


;<")  = s u(n)  ♦ k 


must  be  computed  as  part  of  the  next  SSOR-SI  iteration  u (see  (3)),  the 

pseudo-residual  vector  is  essentially  obtained  as  a byproduct  of  applying  the 
SSOR-SI  method. 

V.  ADAPTIVE  METHOD 

By  using  Theorems  1,  2,  and  4 as  a foundation  we  are  able  to  present 
the  basic  structure  for  an  adaptive  procedure.  For  detailed  descriptions 
of  this  method  and  several  possible  variations  for  the  adaptive  accelera- 
tion of  SSOR,  see  Benokraitis  (1974). 9 

We  state  the  steps  of  the  "algorithm"  in  outline  form  with  a synopsis  of 
the  controlling  theorem(s)  in  the  heading.  We  use  the  word  "algorithm"  loosely, 
since  admittedly  much  is  left  unspecified. 

I.  Theorem  2.  For  any  v ^ 0,  4>(a),v)  < S(S  ). 

— (0 

1.  Choose,  v^,  v2  0. 


2.  Observe 


a.  ^ = (Ku.Vj)  <_  S(Sw) 


b.  d>2  = <J>(u>,v2)  <_  SfS^) 


c.  ^(w)  = max  (4^,  <J>2)  <_  S (S  ) 

VV2 


3.  Minimize  4^(o))  with  respect  to  to  to  obtain  estimate  ojj, 


— — - — - - - — 


. i.-  • — «,  ..'.i3 
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4.  Choose  4,(0),)  as  S (S  ),  an  estimate  of  S(S  ),  (The  situation 
is  depicted  in  Figure  1.) 

II*  Theorem  4.  The  pseudo-residual  vector  6^^  converges  in  direction  to 
dominant  eigenvector  v. 


1.  Set  i = 1 

2.  Iterate  n times  with  SSOR-SI  with  parameters  o).,  S (S  )•  test 

l E to. 

for  convergence  1 

3.  Compute 

6(n)  = n(n)  + k } _ u(n) 

U)  0) 

which  approaches  dominant  eigenvector. 


III. 


4.  Check  if  parameters  should  be  changed. 

a.  If  $(u>j,6^  ^)  < Sg  (S^  ) do  not  change  parameters.  Go  to  1 1.2 

i 

b.  If  $(uk,6^  ^)  > Sg  (S  ) continue  to  step  III  to  change  parameters. 

i 

Theorems  1,  2,  4.  As  6*-11-*  approaches  dominant  eigenvector,  tj>(<*).,  6^) 
approaches  S (S  ). 


1.  Set  vi+2  = 6(n:) 


^i+2  = W 

2.  Observe  4.+1(u)  = max  (4,,  a ....  a ) 

vk,k=l,  ...,  i+2  1 2 

<S(V 

3.  Minimize  4.  . (to)  with  respect  to  to  to  obtain  next  estimate  to. 

11  l + l 

4.  Choose  4.+1(toi+1)  as  S£  (S  ),  an  estimate  of  S (S  ).  Set  i = i+1. 

i+1  “i+1 

Go  to  II. 2.  Process  is  continued  until  convergence. 
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We  briefly  discuss  how  to  choose  n in  step  I I. 2.  Here  we  make  use  of 
the  average  and  asymptotic  average  rates  of  convergence  for  the  SSOR-SI 
method  (Young  (1971)). **  A strategy  which  produces  acceptable  results  is  to 
choose  n so  that  the  average  rate  of  convergence  after  n iterations  is  90%  of 
the  asymptotic  average  rate  of  convergence.  The  convergence  rates  are  computed 
using  the  latest  estimate  of  SfS^).  That  is,  n is  chosen  to  be  the  least  n 

which  satisfies 


• ?loe72ls>  -9  (-108  r) 

1+r 


where 


1 " 


1 + /1-S_(S  ) 
Ev  a) 


A word  about  the  additional  work  required  in  the  adaptive  algorithm  is 
in  order.  Mainly,  the  added  expense  comes  in  changing  the  parameters.  This 
involves  the  computation  of  a and  8,  two  quotients  of  inner  products  in  the 
formula  for  <J> (u> , v)  given  in  (4) -(5).  For  problems  of  the  type  discussed  in 

2 

Section  VI, to  compute  a and  8 requires  approximately  28J  arithmetic  operations 

2 

if  the  mesh  size  is  h = 1/J.  One  SSOR-SI  iteration  requires  approximately  39J 
operations.  Therefore,  four  parameter  changes  are  approximately  equivalent  to 
three  SSOR-SI  iterations  in  terms  of  work  performed.  Since  no  more  than  four 
parameter  changes  were  required  for  the  problems  considered,  the  number  of 
iterations  for  the  adaptive  algorithm  should  effectively  be  increased  by  about 
three  iterations. 


VI.  NUMERICAL  EXAMPLES 

We  present  results  for  a sample  of  the  generalized  Dirichlet  prob- 
lems considered.  The  results  are  given  in  graphic  form  in  Figures  2-4. 

In  each  figure,  we  give  the  differential  equation,  the  region  considered 
and  the  boundary  values.  We  replace  the  differential  equation  by  a 5- 
point  symmetric  difference  eouation  (see  Young  (1977)). 6 

The  number  of  iterations  required  for  varying  mesh  sizes  is  re- 
corded for  optimum,  adaptive,  and  estimated  SSOR-SI  parameters.  In  the 
adaptive  case,  the  subscript  on  the  number  of  iterations  indicates  the 
number  of  parameter  changes  required.  The  estimated  parameters  are  the 
values  of  Young  (1977)6  which  depend  on  bounds  for  the  eigenvalues  of  B and 
LU.  (For  Problem  3.  the  results  for  the  estimated  Darameters  are  not 
given  since  an  excessive  number  of  iterations  are  required.)  The  slopes  s 
of  the  lines  indicate  that  the  number  of  iterations  required  for  convergence 

increases  like  h S. 
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For  smooth  and  some  discontinuous  coefficients  (Problems  1 and  2),  the 

number  of  iterations  required  behaves  like  h an  order-of-magnitude  better 
than  SOR  or  SSOR.  For  cases  involving  higher  discontinuity  (Problem  3),  the 
-3/4 

behavior  is  like  h 
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